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Abstract. We propose a population model for TB-HIV/AIDS coinfection 
transmission dynamics, which considers antiretroviral therapy for HIV infection 
and treatments for latent and active tuberculosis. The HIV-only and TB-only 
sub-models are analyzed separately, as well as the TB-HIV/AIDS full model. 
The respective basic reproduction numbers are computed, equilibria and sta¬ 
bility are studied. Optimal control theory is applied to the TB-HIV/AIDS 
model and optimal treatment strategies for co-infected individuals with HIV 
and TB are derived. Numerical simulations to the optimal control problem 
show that non intuitive measures can lead to the reduction of the number of 
individuals with active TB and AIDS. 


1. Introduction. According with the World Health Organization (WHO), the hu¬ 
man immunodeficiency virus (HIV) and mycobacterium tuberculosis are the first 
and second cause of death from a single infectious agent, respectively [48]. Ac¬ 
quired immunodeficiency syndrome (AIDS) is a disease of the human immune sys¬ 
tem caused by infection with HIV. HIV is transmitted primarily via unprotected 
sexual intercourse, contaminated blood transfusions, hypodermic needles, and from 
mother to child during pregnancy, delivery, or breastfeeding [37]. There is no cure 
or vaccine to AIDS. However, antiretroviral (ART) treatment improves health, pro¬ 
longs life, and substantially reduces the risk of HIV transmission. In both high- 
income and low-income countries, the life expectancy of patients infected with HIV 
who have access to ART is now measured in decades, and might approach that of 
uninfected populations in patients who receive an optimum treatment (see [12] and 
references cited therein). However, ART treatment still presents substantial limi¬ 
tations: does not fully restore health; treatment is associated with side effects; the 
medications are expensive; and is not curative. Following UNAIDS global report 
on AIDS epidemic 2013 [45], globally, an estimated 35.3 million people were living 
with HIV in 2012. An increase from previous years, as more people are receiving 
ART. There were approximately 2.3 million new HIV infections globally, showing 
a 33% decline in the number of new infections with respect to 2001. At the same 
time, the number of AIDS deaths is also declining with around 1.6 million AIDS 
deaths in 2012, down from about 2.3 million in 2005. 

Mycobacterium tuberculosis is the cause of most occurrences of tuberculosis (TB) 
and is usually acquired via airborne infection from someone who has active TB. 
It typically affects the lungs (pulmonary TB) but can affect other sites as well 
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(extrapulmonary TB). In 2012, approximately 8.6 million people fell ill with TB 
and 1.3 million people died from TB. Nevertheless, TB death rate dropped 45 per 
cent between 1990 and 2012, and 22 million lives were saved through use of strategies 
recommended by WHO [48]. 

Individuals infected with HIV are more likely to develop TB disease because 
of their immunodeficiency, and HIV infection is the most powerful risk factor for 
progression from TB infection to disease [18]. In 2012, 1.1 million of 8.6 million 
people who developed TB worldwide were HIV-positive. The number of people 
dying from HIV-associated to TB has been falling since 2003. However, there were 
still 320 000 deaths from HIV-associated to TB in 2012, and further efforts are 
needed to reduce this burden [48] . ART is a critical intervention for reducing the risk 
of TB morbidity and mortality among people living with HIV and, when combined 
with TB preventive therapy, it can have a significant impact on TB prevention [48] . 

Collaborative TB/HIV activities (including HIV testing, ART therapy and TB 
preventive measures) are crucial for the reduction of TB-HIV coinfected individuals. 
WHO estimates that these collaborative activities prevented 1.3 million people from 
dying, from 2005 to 2012. However, significant challenges remain: the reduction of 
tuberculosis related deaths among people living with HIV has slowed in recent 
years; the ART therapy is not being delivered to TB-HIV coinfected patients in the 
majority of the countries with the largest number of TB/HIV patients; the pace of 
treatment scale-up for TB/HIV patients has slowed; less than half of notified TB 
patients were tested for HIV in 2012; and only a small fraction of TB/HIV infected 
individuals received TB preventive therapy [45]. 

The study of the joint dynamics of TB and HIV present formidable mathematical 
challenges due to the fact that the models of transmission are quite distinct [36]. 
Some mathematical models have been proposed for TB-HIV coinfection (see, e.g., 
[2, 3, 22, 28, 30, 36, 40]). In this paper, we propose a new population model for 
TB-HIV/AIDS coinfection transmission dynamics, where TB, HIV and TB-HIV 
infected individuals have access to respective disease treatment, and single HIV- 
infected and TB-HIV co-infected individuals under HIV and TB/HIV treatment, 
respectively, stay in a chronic stage of the HIV infection. 

Optimal control is a branch of mathematics developed to find optimal ways to 
control a dynamic system [10, 16, 31]. While the usefulness of optimal control 
theory in epidemiology is nowadays well recognized (see, e.g., [4, 26, 27, 33, 34]), 
and has been applied to TB models (see, e.g., [5, 15, 19, 21, 41, 42]) and HIV models 
(see, e.g., [23, 29]), to our knowledge optimal control have never been applied to a 
TB-HIV/AIDS coinfection model. In this paper, we apply optimal control theory 
to our TB-HIV/AIDS model and study optimal strategies for the minimization of 
the number of individuals with TB and AIDS active diseases, taking into account 
the costs associated to the proposed control measures. 

The paper is organized as follows. The model is formulated in Section 2. The 
HIV-only and TB-only sub-models of the full TB-HIV/AIDS model are analyzed in 
Section 3 and the full TB-HIV/AIDS model is analyzed in Section 4. In Section 5 we 
propose an optimal control problem and apply the Pontryagin maximum principle to 
derive its solution. In Section 6 numerical simulations and discussion of the results 
are carried out for the optimal control problem associated to the TB-HIV/AIDS 
model. We end mentioning some possible future work in Section 7. 
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2. Model formulation and basic properties. The model subdivides the human 
population into eleven mutually-exclusive compartments, namely susceptible indi¬ 
viduals (S), TB-latently infected individuals, who have no symptoms of TB disease 
and are not infectious (Lt), TB-infected individuals, who have active TB disease 
and are infectious (It), TB-recovered individuals (R), HIV-infected individuals with 
no clinical symptoms of AIDS (Ih), HIV-infected individuals under treatment for 
HIV infection ( Ch ), HIV-infected individuals with AIDS clinical symptoms (A), 
TB-latent individuals co-infected with HIV (pre-AIDS) (Lth), HIV-infected indi¬ 
viduals (pre-AIDS) co-infected with active TB disease (Ith), TB-recovered individ¬ 
uals with HIV-infection without AIDS symptoms (Rh), HIV-infected individuals 
with AIDS symptoms co-infected with active TB (At). The total population at 
time t, denoted by N(t), is given by 

N(t) = S(t) + Lt(I) + It{P) T R(t) + Ih(L) V Cn(t) + A(t) 

+ Ith(L) + L T H(t) + Rh(L) + A T (t). 

The susceptible population is increased by the recruitment of individuals (assumed 
susceptible) into the population, at a rate A. All individuals suffer from natural 
death, at a constant rate /i. Susceptible individuals acquire TB infection from 
individuals with active TB at a rate A t, given by 

'Vr(t) = (Irtt) + lTH(t) + A T (t)), (1) 

where Bi is the effective contact rate for TB infection. Similarly, susceptible indi¬ 
viduals acquire HIV infection, following effective contact with people infected with 
HIV at a rate A h, given by 

Bn 

A H (t) = y [/^(i) + Ith( t) + L T H(t) + Rn(t) + r]c Cn(t) + (A(t) + A T (t ))], 

( 2 ) 

where Bi is the effective contact rate for HIV transmission. The modification pa¬ 
rameter 17^4 > 1 accounts for the relative infectiousness of individuals with AIDS 
symptoms, in comparison to those infected with HIV with no AIDS symptoms. In¬ 
dividuals with AIDS symptoms are more infectious than HIV-infected individuals 
(pre-AIDS) because they have a higher viral load and there is a positive correlation 
between viral load and infectiousness [11]. On the other hand, r\c < 1 translates 
the partial restoration of immune function of individuals with HIV infection that 
use correctly ART [12]. 

Remark 1. For the basic and classical SIR model, one has the force of infection F 
given by F = f3I , which models the transition rate from the compartment of suscep¬ 
tible individuals S to the compartment of infectious individuals I. However, for large 
classes of communicable diseases, it is more realistic to consider a force of infection F 
that does not depend on the absolute number of infectious, but on their fraction with 
respect to the total population N, that is, F = fj jj ■ In our case, the force of infection 
for the HIV is given by X H = ^ [Ih + Ith + Lth + Rh +VcC H + (A + A T )\. 

Only approximately 10% of people infected with mycobacterium tuberculosis de¬ 
velop active TB disease. Therefore, approximately 90% of people infected remain 
latent. Latent infected TB people are asymptomatic and do not transmit TB [43]. 
Individuals leave the latent-TB class Lt by becoming infectious, at a rate fci, or 
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recovered, with a treatment rate T\. The treatment rate for active TB-infected indi¬ 
viduals is t2 . We assume that TB-recovered individuals R acquire partial immunity 
and the transmission rate for this class is given by /3[\t with /3[ < 1. Individ¬ 
uals with active TB disease suffer induced death at a rate dr- We assume that 
individuals in the class R are susceptible to HIV infection at a rate A h- On the 
other hand, TB-active infected individuals It are susceptible to HIV infection, at a 
rate 5Xh, where the modification parameter 5 > 1 accounts for higher probability 
of individuals in class It to become HIV-positive. HIV-infected individuals (with 
no AIDS symptoms) progress to the AIDS class A at a rate p\, and to the class 
of individuals with HIV infection under treatment Ch at a rate </>. Individuals 
in the class Ch leave to the class Ih at a rate wi. HIV-infected individuals with 
AIDS symptoms are treated for HIV at the rate or and suffer induced death at a 
rate d A - Individuals in the class Ih are susceptible to TB infection at a rate i/’At, 
where ^ > 1 is a modification parameter traducing the fact that HIV infection is 
a driver of TB epidemic [24]. HIV-infected individuals (pre-AIDS) co-infected with 
TB-disease, in the active stage Ith , leave this class at a rate p2 ■ A fraction p of 
Ith individuals take simultaneously TB and HIV treatment and a fraction q of Ith 
individuals take only TB treatment. Individuals in the class Ith progress to the 
class Ch at a rate pp2 and to the class Rh at a rate qp2 ■ Individuals in the class 
Ith that do not take any of the TB or HIV treatments progress to the class At at 
a rate (1 — (p + q))p2 , and suffer TB induced death rate at a rate c?t- Individuals 
leave Lth class at a rate 73 . A fraction r of Lth individuals take simultaneously 
TB and HIV treatment and a fraction 1 — r take only TB treatment. Individuals 
in the class Lth progress to the class Ch at a rate rr 3 and to the class Rh at a 
rate (1 — r)r 3 . Individuals in the class Lth are more likely to progress to active TB 
disease than individuals infected only with latent TB. In our model, this progression 
rate is given by £; 2 . Similarly, HIV infection makes individuals more susceptible to 
TB reinfection when compared with non HIV-positive patients. The modification 
parameter associated to the TB reinfection rate, for individuals in the class Rh, is 
given by /3 2 , where /3 2 > 1. Individuals in this class progress to class A, at a rate w 2 . 
HIV-infected individuals (with AIDS symptoms), co-infected with TB, are treated 
for HIV, at a rate 0.2- Individuals in the class At suffer from AIDS-TB coinfection 
induced death rate, at a rate dj^A- The aforementioned assumptions result in the 
system of differential equations 

'S(f) = A - A T (t)S(t) - A H {t)S(t) - pS(t), 

L T (t) = A T (t)S(t) + / 3 [X T (t)R{t ) - (fci + ti + p)L T {t), 

ir(t) = Lt{R) — (t 2 + dT + p + SXnit)) ir(t), 

R(t) = ti L T (t) + r 2 / T (f) - (/?iA T {t) + X H it) + p)R(t), 

i H {t) = A H(t)S(t) — (pi + 4> + ipXT(t) + p)In{t) + a\A(t) + XH{t)R(t) + wiCA(f), 

< A(t) = p\I H {t) + (t) - A(t) - (p + d A )A(t), 

Ch (t) = 4*1 h (t) + PP2lTn(t) + rT 3 L T H(t ) - (wi + p)Cn{t), 

Lth^) = f3 2 X T (t)RH{t ) — ( k 2 + t 3 + p)L T H{t), 

irH{t) = SXH(t)lT(t) + ^Xt( t)I h( t) + a 2 Ar(t) + k 2 LTH{t ) — (P 2 + P + d T ) Ith{$), 

Rii(t) = qp2lTn{t) + (1 — r) r 3 LTH{t) — ^/ 3 2 A-r(t) + w 2 + pj Rni^t), 

At{ t) = (1 — {p + q))p2lTH{t) — (a 2 + P + dTA)Ar(t), 


( 3 ) 
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Figure 1. Model for TB-HIV/AIDS transmission. 


that describes the transmission dynamics of TB and HIV/AIDS disease. The model 
flow is illustrated in Figure 1. 


2.1. Positivity and boundedness of solutions. Since the system of equations 
(3) represents human populations, all parameters in the model are non-negative and 
it can be shown that, given non-negative initial values, the solutions of the system 
are non-negative. Consider the biologically feasible region 

^ = {{S,L t ,It,R^Ih,A,Ch,Lth,Ith,Rh,At)€:^ 1 ^ : N< A/p}. 

In what follows we prove the positive invariance of (i.e., all solutions in Q remain 
in Q for all time). The rate of change of the total population, obtained by adding 
all the equations in model (3), is given by 


dN 

dt 


A — fiN(t ) — dT^r(i) — dAA(t) — dTlTuif) ~ dTAAx(t). 


Using a standard comparison theorem [25] we can show that 


N(t) < N(0)e~^ 


+ - (1 - e~^) . 


In particular, N(t) < ]”) if iV(0) < —. Thus, the region is positively invariant. 
Hence, it is sufficient to consider the dynamics of the flow generated by (3) in f2. 
In this region, the model is epidemiologically and mathematically well posed [20]. 
Thus, every solution of the model (3) with initial conditions in Q remains in Q for 
all t > 0. This result is summarized below. 


Lemma 2.1. The region H is positively invariant for the model (3) with non¬ 
negative initial conditions in Ml 1 . 
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3. Analysis of the sub-models. In this section we analyze the models for HIV 
only (HIV-only model ) and TB only (TB -only model). 

3.1. HIV-only model. The model that considers only HIV (obtained by setting 
Lt = It = R = Lth = Ith = Rh = At) is given by 

’S(t) = A-\ H (t)S(t)-pS(t), 

lH{t) — A u{t)S(t ) — (pi + 4> + ii)Iu(t) + ai A(t) + u>iCh, 

< . ( 4 ) 

A(t) = PiIh(1) — ( oi\ + P + d A )A(t), 

— (wi + p)Cii{t), 

where 

A H(t) = [IHit) + r)cC H (t) + VAA(t)] 

with 

N(t) = S(t) + I H {t) + A(t) + C H (t). 

Analogously to Lemma 2.1, we can prove that the region 

ni = {iS 7 I H ,A,C H )eRl : N<A/p} (5) 

is positively invariant and attracting. Thus, the dynamics of the HIV-only model 
will be considered in Hi. 

3.1.1. Persistence. In this section, we look for the conditions under which the host 
population and disease will persist. Rewriting the submodel system (4) as 

'S[t) = A - PAN)(Ih+vcC h+va a) s ^ _ 

i H (t) = dAN)(i H +vcC H+VA A) s{t) ^ (pi + ^ + ^ (t) + ai A (t) + uiCh, 

< 

Ch it) = <1)1 h it) — (wi + p)Cnit), 

.A(t) = pJHit) - iai+p + d A )A(t), 

( 6 ) 

in what follows we assume that /3 2 (AT) is continuous for N > 0 and continuously 
differentiable for N > 0; /3 2 (A) is monotone nondecreasing in N; and /3 2 (1V) > 0 if 
N > 0. 

Remark 2. In this work /3 2 denotes the effective contact rate for HIV transmission. 
It is a constant for a concrete situation, but one can look to it as variable in the 
sense that, depending on the situation/region, one can have different values for this 
parameter. This is so because /3 2 is related with the level of contagion/propagation 
of the disease. In Section 5 we consider fixed values for fdi and /? 2 , which represent 
specific cases of the infection level. By varying j3i and /? 2 we vary the basic repro¬ 
duction numbers (see expressions (12) and (20) for Ri and 1? 2 , respectively). Here 
we consider /? 2 as a function of N to discuss persistence. 

It is convenient to reformulate the model in terms of the fractions of the S, Ih , 
A and Ch parts of the population, 

S Ih C h A 
x =y = ~Vr > z =^n W =S~r 


(7) 
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and express (6) in these terms to obtain the system 
'N = A - (p + d A w)N, 

x= ^(1-x) - /3 2 (N)(y + r) C z + r) A w)x + - y + d A w), 

< y = /3 2 (N)(y + r] C Z + y A w)x + yd A w - (pi + </> + y) y + ai w + uqz, (8) 
z = <fy — (w i + y — d A w)z, 

,ai = piy - («i + y + d A - d A w)w. 

Equations (7) suggest that x + y + z + w = 1. The manifold x + y + z + w = 1, 
x , y,z,w > 0, is forward invariant under the solution flow of (8), which has a global 
solution satisfying (7). We now show conditions under which the host population 
will persist. 

Theorem 3.1. Let (3 2 ( 0) = 0, iV(0) > 0. Then the population is uniformly persis¬ 
tent, that is, 

liminf N(t) > e, 

t—¥ OO 

where e > 0 does not depend on the initial data. 

Proof. We have to show that the set 

X 2 = {N = 0, x>0,y>0, z>0,w>0, x + y + z + w = l} 
is uniform strong repeller for 

X\ = {N >0, x>0,y>0,z>0, w>0,x + y + z + w=l}. 

Theorem 3.2, Theorem 3.3 and Corollary 1 are taken from [3, 44], 

Theorem 3.2. Let X be a locally compact metric space with metric d. Let X be 
the disjoint union of two sets Xi and X 2 such that X 2 is compact. Let be a 
continuous semiflow on X\. Then X 2 is a uniform strong repeller for X\, whenever 
it is a uniform weak repeller for Xi. 

Theorem 3.3. Let D be a bounded interval in K and g : (to,oo) x D —> R be 
bounded and uniformly continuous. Further, let x : ( t 0 ,oo ) —> D be a solution of 

x' = g{t,x), 

which is defined on the whole interval (to,oo). Then there exist sequences s n ,t n —> 
oo such that 

lirn 5 ( Sn ,Xoo)=0 = lim g(t n ,x°°). 

n—f oo n—foo 

Corollary 1. Let the assumptions of Theorem 3.3 be satisfied. Then 

a) lim inf^oo g(t, x^) > 0 > limsup^^ g{t, x^), 

b) lim inft-^oo g(t, x°°) > 0 > lim sup g(t,x°°). 

As the assumptions of Theorem 3.2 are satisfied, it suffices to show that X 2 is a 
uniform weak repeller for X\. Let r = y + z w. Then, 

r' = f 2 {N)(yx + rjczx + rj A wx) + d A wr — pr — d A w 

- + d A (r - 1), 


< (32(N)(1 + r)c + Va) 
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using the fact that x, y , z,w,r < 1. This implies that 

_A_ r oo + (1 _ r ~ )d A < /3 2 (JV°°)(1 + ^ + VA ) 


fciN 00 ) > 


A r° 


(1 — r°°)d A 


( 9 ) 


TV 00 (1 + vc + rj A ) 1 + Vc + Va' 
From the equation of N in ( 8 ) we have 

liminf > - 7 ^-(/i + d A w°°) > —^- (y + d A r °°). 

t->oo N dt - N°° - N°° ^ ' 

Hence N increases exponentially, unless 

A . 1 / A 


N°° 


< y + d A r°°, that is, 


d A \N 


-y)<r° 


( 10 ) 


Combining (9) and (10), we obtain that 
h(N°°) > ' A 


1 


d A N°°(l +r] C + Va) 1 + ilc + VA 


A 

N°° 




d A 


1 + Vc + Va 

( 11 ) 

as /? 2 (0) = 0 and /3 2 (iV) is continuous at 0, A ' 00 > £ > 0 with e not depending on 
the initial data. From (11) we see that we can relax /3 2 (0) = 0 and require 


182(0) < 


A 


^d A N°°(l + r)c + Va) 1 +t?o+i?a 
This concludes the proof. 


A 

N°c 


d A 


1 + yc + va 


□ 


The disease is persistent in the population if the fraction of the infected and 
AIDS cases is bounded away from zero. If the population dies out and the fraction 
of the infected and AIDS remains bounded away from zero, we would still say that 
the disease is persistent in the population. 

Proposition 3.4. Let fa (00) (I + tjc + Va) > j^r°°. Then the disease is uniformly 
weakly persistent insofar as 

r°° = limsupr(t) > e, 

t—t OO 

with e > 0 being independent of the initial data, provided that r(0) > 0. 

The proof of Proposition 3.4 is outlined in [44]. 

3.1.2. Local stability of disease-free equilibrium. The model (4) has a disease-free 
equilibrium (DFE), obtained by setting the right-hand sides of the equations in the 
model to zero, given by 

z 0 = (s*,r H ,A*,cr H ) = 0 , 0 , 0 ). 

The linear stability of E 0 can be established using the next-generation operator 
method on the system (4) . Following [46] , the basic reproduction number is obtained 
as the spectral radius of the matrix FV^ 1 at the DFE, So, with F and V given by, 
respectively, 

0 0 0 

\ TT fa s fa a a 8 

N N 

0 0 0 

0 0 0 


0 

02 VC s 

N 

0 

0 


F = 
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and 

r \ | 02 s 02VaS 02VC s 1 

P ~JT N N 

0 Cl — Q. i — 

v= , 

0 -p 1 c 2 o 

0 —<j) 0 c 3 

where Ci = p\ + cf + p, C 2 = ay + p + cLa, C 3 = u>i + p. The basic reproduction 
number is given by the dominant eigenvalue of the matrix FV -1 , that is, 


fa A [C 3 (C 2 + r]A pi) + pc <t>C 2 ) 

Np [p (C 3 (pi + C 2 ) + C 2 cf> + Pi^a) + Pi^icIaI 


( 12 ) 


The basic reproduction number R\ represents the expected average number of new 
HIV infections produced by a single HIV-infected individual when in contact with 
a completely susceptible population [46]. 


Remark 3. The next-generation matrix is one of the most well known methods 
in epidemiology to compute the basic reproduction number for a compartmental 
model of the spread of infectious diseases. To calculate the basic reproduction 
number through this method, the whole population is divided into n compartments 
in which there are m < n infected compartments. Let Xi, i = 1,2,..., to, be the 
numbers of infected individuals in the ith infected compartment at time t. Now, 
the epidemic model is x[ = F.i(x) — Vi(x) or, in vector form, x' = F{ x) — V(x). Let 
Xq denote here the disease-free equilibrium state. The Jacobian matrices of F(x) 
and V(x) are, respectively, 


DF{x) 


F 0 
0 0 


and DV (x) 


V 0 
fa fa ’ 


where F and V are the m x m matrices given by 


dFj(x q) 
dxj 


dVi(x Q ) 

dxj 


The matrix FV^ 1 is known as the next-generation matrix and its spectral radius is 
the basic reproduction number of the model. The reader interested in all the details 
about the computation of the basic reproduction number by the next-generation 
matrix is referred to [14, 46] or any good book on dynamical modeling and analysis 
of epidemics (e.g., [13]). 


Lemma 3.5. The disease free equilibrium So is locally asymptotically stable if Ri < 
1, and unstable if R\ > 1. 


Proof. Following Theorem 2 of [46], the disease-free equilibrium, So, is locally 
asymptotically stable if all the eigenvalues of the Jacobian matrix of the system 
(4), here denoted by M (So), computed at the DFE So, have negative real parts. 
The Jacobian matrix of the system (4) at disease free equilibrium S 0 is given by 


M (S 0 ) 


—M 


02 A 

fiN 


02 UA A 02 Vc A 

fi N fiN 


0 


02 A _ si 
UN 


02 VA A 

UN 


+ Oil 


02 VC A 

fj, N 


+ Wi 


0 pi c 2 


0 


0 


C 3 


(13) 


0 
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One has 

trace [M (So)] = — P + ^ ^ ~~ (Ci + 6*2 + 6*3) < 0 

and 

det [M (Sq)] = — jj- [ 03 ( 6*2 + PiVa) + (fiVcC 2 ] 

+ /x [/x (6*3(pi + 6*2) + C*20 + picIa) + > 0 

for i?i < 1. We have just proved that the disease free equilibrium S 0 of model (4) 
is locally asymptotically stable if R\ < 1, and unstable if R\ > 1. □ 


3.1.3. Global stability of disease-free equilibrium (DFE). Following [8], let us rewrite 
the submodel system (4) as 


dX 

dt. 

dZ 

dt 


= F(X, Z), 
= G(X,Z), 


G(X, 0) = 0, 


(14) 


where X = S and Z = ( Ih , A, Ch), with X £ R + denoting the total number of un¬ 
infected individuals and Z £ denoting the total number of infected individuals. 
The disease-free equilibrium is now denoted by 


U o = (X o ,0), 


where X 0 = 



The conditions (HI) and (H2) below must be met to guarantee global asymptotically 
stability: 

(HI) for = F(X, 0), Uo is globally asymptotically stable; 

(H2) G(X, Z) = AZ-G[X , Z), G{X, Z) > 0 for (X, Z) £ G, where A = D Z G(U 0 ,0) 
is a Metzler matrix (the off diagonal elements of A are nonnegative) and Q is 
the region where the model makes biological sense. 


Theorem 3.6. The fixed point Uq = (Xo,0) is a globally asymptotically stable 
equilibrium of (4) provided R\ < 1 and the assumptions (HI) and (H2) are satisfied. 


Proof. We have 

^ = F(X, Z) = [ A - X H S — pS ] , 

F(X, 0)= [ A -pS ] , 

A H{t)S(t ) — C 1 I H (t) + aiA(t) + ujiCh 
P iI H (t) - C 2 A(t) 

- C 3 C H (t) 

and G(X, 0) = 0. Therefore, 


dZ 

dt 


= G(X, Z) = 


dX 

dt 


F(X, 0) 


A — pS 
—pRr 


hi A 02 VA A , 02 TIC A , , , 

TvjT 01 - l_ai ■VA-rWl 


N/i 

-C 2 


Nn 


0 -6*3 


A = D Z G(X 0 , 0) 


Pi 
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and 


G(X,Z) = 

' G 3 {X,Z) ' 

G 2 (X,Z) 

— 

P 2(1 - + VaA + tjcCh) 

0 


. G 3 (X,Z) _ 


0 


(15) 


It follows that Gx(X,Z) > 0, G 2 (X,Z) = G 3 {X,Z) = 0. Thus, G(X,Z) > 0. 
Conditions (HI) and (H2) are satisfied, and we conclude that Uq is globally asymp¬ 
totically stable for i?! < 1 . □ 


3.1.4. Existence of an endemic equilibrium. To find conditions for the existence of 
an equilibrium for which HIV is endemic in the population (i.e., at least one of Ifj, 
A* or Cff is non-zero), denoted by E h = (S*, Ifj, A *the equations in (4) are 
solved in terms of the force of infection at steady-state (X* H ), given by 


A 


* _ 

H ~ 


@2 {Ijj + VaA* + rjc Cjj) 
N* 


(16) 


Setting the right hand sides of the model to zero (and noting that A h = A * H at 
equilibrium) gives 


S* = 


A 


A *h + P 


T* — _ 

1 H — 


X* H AC 2 C 3 

D 


A* = -- 


Pi\* h AC 3 


D 


C* H = - 


<t>X* H AC 2 


' H 

D 


(17) 

with D = -(X* H + (C 3 (pi + C 2 ) + C 2 4> + p\dA) + Pi^idA)- Using (17) in the 
expression for A^ in (16) shows that the nonzero (endemic) equilibria of the model 
satisfy 

.* _ A/3 2 (C 2 C 3 + rjAPiC 3 + rjcf’G'z) 

H Al [p (C 3 (pi + C 2 ) + C 2 <j> + Pi 6 ?a) + Pi^hcUi] 


that is, 


X* H = -p(l-R 1 ). 

The force of infection at the steady-state X* H is positive, only if Ri > 1. We have 
just proved the following result. 


Lemma 3.7. The submodel system (4) has a unique endemic equilibrium whenever 
Ri > 1. 


3.1.5. Local stability of the endemic equilibrium. In what follows we prove the local 
asymptotic stability of the endemic equilibrium E h, using the center manifold the¬ 
ory [6], as described in [9, Theorem 4.1], with E h = {S*,Ifj, A*,Cjj) and each of its 
components given as in (17). To apply this method, the following simplification and 
change of variables are first made. Let S = aq, Ih = x 2 , A = x 3 and Cjj = X 4 , so 
that N = xi +x 2 +x 3 + X 4 . Further, by using vector notation X = (x 3 , x 2 , x 3 , X 4 ) T , 
the submodel (4) can be written in the form dZL = (/ 1; / 2 , / 3 , / 4 ) r , as follows: 

dx 1 /f 2 [x 2 + rj A x 3 + r} C x 4 ) 

—— = fi = A - Xi — LIX 1, 

dt x\ + x 2 + x 3 + X4 


dx 2 
dt 


fd 2 {x 2 + p A x 3 + pcx 4 ) 
Xi + x 2 + x 3 + £4 


xi — C 1 X 2 + a 3 x 3 + W 1 X 4 , 


dx 3 
dt 


= h = P\x 2 - C 2 x 3 , 


dx4 

dt 


f4 = 4>X 2 — C 3 X4- 


(18) 
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The basic reproduction number of the submodel (4) is given by (12). Choose as 
bifurcation parameter fa, by solving for fa from R\ = 1: 

o* _ P (^3 (Pi + C 2 ) + C 2 (f> + Pld A ) + PlU\d,A 
C 3 {C 2 + pa Pi) + Pc c/>C 2 

The submodel (4) has a disease free equilibrium given by 

S H = (^ 10 ,^ 20 ,^ 30 , 2 : 40 ) = (—,0,0,0 

\P 

The Jacobian of the system (18), evaluated at £ 0 , M (£ 0 ), and with fa = fa, is 
given by (13). Note that the above linearized system, of the transformed system 
(18) with fa = fa, has a zero eigenvalue which is simple. Hence, the center manifold 
theory [6] can be used to analyze the dynamics of (18) near fa = fa. In particular, 
Theorem 4.1 in [9] is used to show the locally asymptotically stability of the endemic 
equilibrium point of (18), for fa near fa. 

The Jacobian M(£ 0 ) at fa = fa has a right eigenvector (associated with the 
zero eigenvalue) given by w = [w\, iv 2 , w 3 , w^\ T , where 

w 3 p(C 3 {p\ + C 2 ) + C 2 <f> -\- p\d a) + piU\d,A 
Wl ~ C 3 PP 1 ’ 

w 3 C 2 
w 2 =-, 

Pi 

w 3 =w 3 >0, 

C 2 (j)'w 3 
WA ~ C 3 p 1 ' 

Further, M(£ 0 ) for fa = fa has a left eigenvector v = [iq, v 2 , v 3 , iq] (associated 
with the zero eigenvalue), where 

Vi = 0, 

C 2 C 3 + pAPlC 3 + 77C0C2 

V<2 — V3 - 

' p{paCi + PAUi + 0 : 1 ) + uj l p A p 1 + ( p C (j) + Wi)«i ’ 

v 3 =v 3 > 0, 

_ PcC 2 + p(uji + pc PI + P c4>) + Pi(ujipa + Pcd A ) + (c 2 i + pcfajai + d A ) 

piPACi+PA^i+a^+ufauPi+a^+pc^ai 

To apply Theorem 4.1 in [9] it is convenient to let fk represent the right-hand side of 
the fcth equation of the system (18) and let Xk be the state variable whose derivative 
is given by the fcth equation for k = 1,..., 4. The local stability near the bifurcation 
point fa = fa is then determined by the signs of two associated constants, denoted 
by a and b, defined (respectively) by 


4 

E 


V k WiWj 


d 2 fk 

dxidx. 


(0,0) and b = 


4 


2- 

k,i= 1 


with (j) = fa — fa. Note that, in f k { 0,0), the first zero corresponds to the DFE, 
S 0 , for the subsystem (4). In other words, fkfafa) = 0, for k = 1,..., 4, if and 
only if the right-hand sides of the equations of (4) is zero at So- Moreover, from 
(j) = fa — fa we have 4 > = 0 when fa = fa, which is the second component in f k { 0, 0). 
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For the system (18), the associated non-zero partial derivatives at the disease 
free equilibrium Eo are given by 

cPj\ = 2(3^i d 2 h _ /3*/r(l + va) g> 2 /i __ p*p(l+p c ) 

dx 2 A ’ dx 2 dx 3 A ’ 8 x 2 8 x 4 A 

d 2 fi _ 2/3* pp A d 2 fi _ /3*p{p A + Vc) 


dx 3 A ’ dx 3 dx 4 A 

d 2 fi _ /3*p(l + Vc ) d 2 h _ 2/3*/i ?? c 


<9x4 <9*2 

A 

’ dx 2 

A 


8 2 h 2/3*/i 

8 2 h 

P*p(i + va) 

8 2 f 2 

/3*p{l + p c ) 

<9x 2 A ’ 

dx 2 dx 3 

A 

8 x 2 8x4 

A 


8 2 h 

2/3*p A p 

8 2 f- 2 

/3*p(va + vc) 

8x\ 

A ’ 

8x 3 8x4 

A 

8 2 f 2 

2/3 *vcp 

8 2 f2 

/3*p(l + vc) 

H 

CO 

A ’ 

<9*49*2 

A 


It follows from the above expressions that 

a _ w 3 C2<l>P 2 H (va + Vc) _ 2 DiC% [C 3 + ^(1 + 77c)) 
pi (wi + p) A D 2 p\C 3 

2Di (C 2 C 3P i (C 3 (l + va) + </>) + Cf»MP? + C 2 ^) 
D 2 p\C 2 

with 

Di = v 3 w 3 (3 2 p ( C 3 C 2 + C 3 t]aPi + Vcffii) > 

D 2 — A {vaPiC 3 + C 3 a 1 + vapC 3 + pc/> + Vc<t>u i). 

For the sign of b , it can be shown that the associated non-vanishing partial deriva¬ 
tives are 

d 2 fi d 2 h d 2 h 

dx 2 d/3* ’ < 9*3 <9/3* ^ Al 8 x 48 / 3 * ^ C ’ 

8 2 h _ d 2 h _ d 2 h 

<9x 2 <9/3* ’ 8x 3 8/3* Va ’ 8 x 48 / 3 * Vc ' 

It also follows from the above expressions that 

k _ ((va + Vc)Pi + C 2 ) (C 3 C 2 + C 3 t]aPi + vefiC 2 ) U 3 IU 3 
{(VApi + ai + vap)C 3 + VAPf> + VC</>ol\) pi 

From the previous computations, we have a < 0 and b > 0. Thus, using Theorem 4.1 
of [9], the following result is established. 

Theorem 3.8. The endemic equilibrium E h is locally asymptotically stable for the 
basic reproduction number R\ (12) near 1. 
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3.2. TB-only model. The sub-model of (3) with no HIV/AIDS disease, that is, 
Ih = A = Ch = L T h = Ith = Rh = A t = 0, is given by 

'S(t) = A - X T (t)S(t) - nS(t), 

L T (t) = X T (t)S(t) + T {t)R(t) - Oi + ri + fi)L T (t), 

. (19) 

It it) = kiLr{t) — (72 + dr + n)Ir(t), 

R(t) = t\L t ( t) + t 2 /t (0 - {P[X T (t) + n)R(t), 


where 


and 


Xr(t) 


thI T {t) 

N{t) 


Nit) — Sit) + Lx(t) + It it) T Rit). 


The sub-model (19) was proposed and analyzed in [7]. This model incorporates the 
basic properties of TB transmission and dynamics. The basic reproduction number 
i?2 of (19) is given by 


A ( pl ) ( kl ] 

UN \ jj, + d T + r 2 / \n + ki + Ti J 


( 20 ) 


The existence, uniqueness and local asymptotic stability of the disease-free and 
endemic equilibria are proven in [7, Theorem 1], 


4. Analysis of the full model. We now consider the full model (3), with the 
DFE given by 


7 _ ( QO T o to do TO AO rto TO TO DO 40 \ 

-o — Ia lj TH , ± th , 

A 


= ( -, 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0,0 


The associated matrices F and V (see Section 3.1.2) are given, respectively, by 


F = [ Fi F 2 \ 

with 


0 

0 

0 

0 

0 

0 

At 

0 

/3i (S+R/3[) 

N 

P'\X T 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

A H 

0 

0 

A H 

0 2(5+1?) 

N 

02Va{S+R) 

N 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

02 01 Rh 

N 

0 

0 

0 

0 

0 

SXh + ^iM 

0 

A^+V’A T 

S02VaIt 

N 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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0 



0 


0 

0 

0 


0 


0i(S+R0() 

N 


0 

N 

0 


0 



0 


0 

0 

0 


0 



0 


0 

0 

02Vc(S+R) 

N 

0 2 (S+R) 

N 

0 2(S+R) 

N 


0 2(S+R) 

N 

02Va(S+R) 

N 

0 


0 



0 


0 

0 

0 


0 



0 


0 

0 

0 

02 01 Rh 

N 

P2 Xt 


02 01 Rh 

N 

0 

S(3 2 It 

N 

S/3 2 It 

N 

+ 

4’0iIh 

N 

5(3 2 It 

N 

80 2 T)aIT 1 iP0iIh 

N ' N 

0 

0 


0 



0 


0 

0 

0 


0 



0 


0 

0 

= [Vi V 2 \ with 








I" Xt + A h + P 0 


PiS 

N 


0 

S 

N 

02 Sr] a 
N 

0 

c 4 


0 



0 

0 

0 

0 

-k\ 


SX H + c 5 


0 

5(3 2 I t 

N 

802V aIt 
N 

0 

-n 


0'i0i R- 
N 

t 2 N 

ft At 

+ A 

H 1 P jy 

02V A R 
N 

0 

0 


it>0i Ih 

N 


0 

'i/jXt H - Ci — cx 1 

0 

0 


0 



0 

-Pi 

c 2 

0 

0 


0 



0 


0 

0 

0 


0 



0 

0 

0 

0 

0 


0 



0 

0 

0 

0 

0 


020lR H 

N 


0 

0 

0 

0 

0 


0 



0 

0 

0 

02VcS 

N 

P 2 5 

N 


(01+0 2)S 

N 


P 2 S 

N 

(01+02 Va)S 

N 

0 

0 



0 



0 

0 

Vc It 

N 

8fi 2 It 
N 



S/3 2 It 

N 



5(3 2 It 

N 

802V aIt 

N 

02VcR 

N 

/ 3 2 R 

N 


(p[pi+p2')R 

N 


0 2 R (0101+0 2Va)R 

N N 

~+>l 

0 



4>0i Ih 

N 


0 

tl>0i Ih 

N 

0 

/3 2 R 

N 



0 



—w 2 

0 

+ fJ, 

-rr 3 



-PP2 


0 

0 

0 

c 6 



0 



0 

0 

0 

-k 2 



c 7 



0 

— 0:2 

0 (- 

1 + r) r 3 

, 0201 Rh 

Q P 2 \ jy 

Z?2^T H - ^2 + fJ> 

P' 2 Pi Rh 

N 

0 

0 


(-1 +p + q) p 2 


0 

c 8 
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where C 4 = k\ + t 4 + p, C 5 = r 2 + p + g?t, + t 3 + Ab C 7 = P 2 + P + dr and 

Cs = 0.2 + cLta + P- The dominant eigenvalues of the matrix FV -1 are 

o __ 02 A (C 3 (C 2 + Pi) + 0^2 ) _ „ _ 

A/ - At [a* {Cs(pi + C 2 ) + C 2 <$> + Pid A ) + Pi^id^] ’ PNC 5 C 4 

Thus, the basic reproduction number i?o of the model (21) is given by 

Ro = max{i?i, R 2 }. 

Using the same procedure as in Section 3.1.2, the following result holds. 

Lemma 4.1. The DFE of the full HIV-TB model (3), given by So, is locally asymp¬ 
totically stable if Rq < 1, and unstable if Rq > 1. 

Remark 4. There are different ways to compute the basic reproduction number 
Rq. Here we are computing it using one of the most well-known methods: Rq 
is the dominant eigenvalue of the associated next-generation matrix FV^ 1 (see 
Remark 3). A justification for the value of the basic reproduction number Rq to be 
max{f?i, i? 2 } is given in Section 4.4 of [46]. 

5. Optimal control problem. In this section we present an optimal control prob¬ 
lem, describing our goal and the restrictions of the epidemic. In the model without 
controls discussed so far, we have p representing the fraction of Ith individuals that 
take HIV and TB treatment and q representing the fraction of Ith individuals that 
take TB treatment only. Roughly speaking, the problem of optimal control consists 
to determine the optimal combination for the values of p and q. For this reason, we 
take p as the control iti and q as the control u 2 . Precisely, we add to the model (3) 
the two control functions ifi(-) and u 2 (-) i n the following way: 

'Sft) = A - A , T {t)Sft) - X H (t)S(t) - a iS{t), 

L T {t) = A T ft) Sft) + fi[X T (t)R{t) - (fci + ti + a«) L T (t), 

1r(t) = kiLr{i) — (r 2 + dx + p + SX n(i)) Irft), 

R(t) = TiL T {t) + r 2 F T (t) - {f)[X T {t) + A H (t) + p)R{t), 

i H (t) = A H (t)S(t) - (pi + (j) + ipX T (t) + p)I H ft ) + aiA(t) + X H (t)R(t) + C H (t), 

< Aft) = pilHit) + w 2 i?jy(t) - aiA(t) - (ai + d A )A(t), 

Cn{t) = ) + u\(t) p2lrH{t) + rT 3 L T H{t) - (wi + p)Cn{t), 

LrHit) = @2 ^t(A)Rh if) — (fe + T 3 + p) L'Tff(t), 

irH(t) = SXH(t)lT(t) + ifX t( t)I h( t) + a 2 H T (t) + fc 2 Lt.fi (t) — (P2 + P + dr) Ith ft), 
Rni^t) = U2ft)P2 ITh ft) + (1 — ) t^Lth if) — (p 2 Xrft) + w 2 + p^j f?n(t), 

At ft) = (1 — (uift) + it 2 (i))) P2.Ith[1) — (02 + P + dr A ) At ft). 

( 2f) 

As already mentioned, the controls u\ and rt 2 represent the fraction of Ith in¬ 
dividuals that are treated for TB and HIV (simultaneously) and treated for TB 
only, respectively. If we consider fixed values for u\ and zi 2 in (21), then we get 
the model (3) with u\ = p and u 2 = q. The aim is to find the optimal values 
u\ and '« 2 of the controls u\ and w 2 , such that the associated state trajectories 
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S*, L?p, 1^1 R*, I*hi A*, Cft , Ltp H , I^ H , R* h , Alp, solution of the system (21) in the 
time interval [0,T] with initial conditions S*(0), Lp( 0), 7£(0), 7?*(0), 7^(0), A*(0), 
0), L?p H ( 0), Irp H (0), R* H (0), A^(0), minimize the objective functional. Here the 
objective functional considers the number of HIV-infected individuals with AIDS 
symptoms co-infected with TB At, and the implementation cost of the strategies 
associated to the controls ttj, i = 1,2. The controls are bounded between 0 and 
0.95. We assume that p and q cannot take values greater than 0.95 because we 
assume that there are some budgetary constraints or some resistance from patients 
in making the treatments (treatment for HIV and TB together or just the treatment 
for TB). In other words, we assume that one cannot treat all the people for both 
diseases or even just for tuberculosis. This is more than reasonable from biological 
side. Moreover, the sum of p + q is also taken as bounded by 0.95. This is related 
with the formulation of the model. Indeed, note that 1 — (j> + q) is the fraction of 
Ith individuals who are not treated for TB and HIV simultaneously and are also 
not treated for TB alone. For this reason, what we assume is that this fraction 
of individuals takes at least the value of 5%. This is in agreement with available 
medical data. Precisely, we consider the state system (21) of ordinary differential 
equations in R 11 with the set of admissible control functions given by 


0 = M-), U 2 (0) e (A°°(0, T)) 2 I 0 < ui(i), U 2 (t) < 0.95 


and 0 < ui(t) + u 2 {t) < 0.95, Vi € [0, T] >. (22) 


The objective functional is given by 

r T 


J{ui{-),u 2 {-)) = [ 
Jo 


A / \ ^ r .1 O / \ IVo O / \ 

A r(t) + -rrUi(i) + -tt u 2 (t ) 


dt , 


(23) 


2 2 

where the constants Wi and W 2 are a measure of the relative cost of the interventions 
associated to the controls ui and u 2 , respectively. 


Remark 5. Epidemiologically, our cost functional tells us that we want to mini¬ 
mize the number of HIV-infected individuals with AIDS symptoms co-infected with 
active TB. For that, one applies control measures that are associated with some 
implementation costs that we also intend to minimize. Other cost functionals may 
be used as well. Here, by considering the cost with controls in a quadratic form, 
we are being consistent with previous works in the literature (see, e.g., [35, 42]). 
Moreover, a quadratic structure in the control has mathematical advantages: if the 
control set is a compact and convex polyhedron (as it is the case here), it imply 
that the Hamiltonian attains its minimum over the control set at a unique point. 
For future work we plan to compare the results now obtained, for a cost with a 
quadratic form in the controls, with those of a linear cost in the controls. 


In order to simplify the formulation of the optimal control problem, let fi rep¬ 
resent the right-hand side of the ith equation of system (21), Xi be the state vari¬ 
able whose derivative is given by the ith component of F, F = (/i,..., /n), and 
X = (aq,..., ecu), i = 1,..., 11. We consider the optimal control problem of deter¬ 
mining X*(-) associated to an admissible control pair (u*(-), u 2 (-)) £ 0 on the time 
interval [0,T], satisfying (21), the initial conditions A(0) and minimizing the cost 
function (23), that is, 

J(ul(-),u 2 {-)) = nun J(ui(-),u 2 (-))- 


(24) 
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The existence of optimal controls (u^(-), u 2 (')) comes from the convexity of the cost 
functional (23) with respect to the controls and the regularity of the system (21) (see, 
e.g., [10, 16] for existence results of optimal solutions). According to the Pontryagin 
maximum principle [31], if («,][(•), u 2 (-)) € 0 is optimal for the problem (21), (24) 
with the initial conditions X(0) and fixed final time T, then there exists a nontrivial 
absolutely continuous mapping A : [0, T] -A R 11 , A(i) = (Ai(t),..., An(i)), called 
adjoint vector , such that 

dH dH 

Xi =—{X,X,u 1 ,u 2 ), Xi =- — (X,X,ui,u 2 ), i = 1, - - -, 11, (25) 

where the function H = H(X, X, ui,u 2 ) defined by 

Wl n Wo a 

H = At H——% ^—2~ U2 u i-> u 2 )) 

is called the Hamiltonian , and the minimality condition 

H{X*(t),X*(t),u* 1 {t),u* 2 (t))= min H(X*(t), A*(i), u u u 2 ) (26) 

0<ii 1 ,u 2 <0.95 

ni+n2<0.95 

holds almost everywhere on [0, T]. Moreover, the transversality conditions 

Xi(T) = 0, i = 1,.. -, 11, 

are also satisfied. 

6. Numerical results and discussion. In this section we present results of the 
numerical implementation of extremal control strategies for the TB-HIV model (21). 
First we solve numerically the optimal control problem (21), (24) with initial condi- 


S( 0) 

66IV(0) 

120 

L t ( 0) 

37AT(0) 

120 

/t(0) 

5JV(0) 

120 

m 

2N(0) 

120 

Mo) 

2N(0) 

120 

A(0) 

N( 0) 
120 

C H { 0) 

Lth{ 0) 

Ith{ 0) 

Rh{ o) 

At( 0) 


N{ 0) 

2AT(0) 

2N(0) 

N( 0) 

N(0) 


120 

120 

120 

120 

120 



Table 1. Initial conditions of the TB-HIV/AIDS model, where 
V(0) = 30000. 


tions given in Table 1, and fixed final time T = 50 years. The initial conditions were 
estimated as follows. We assume that more than half of population (55%) belongs 
to the subgroup of susceptible and that a big percentage (~ 31%) is infected with 
TB but is in the latent stage. This is justified from the fact that “about one-third of 
the world’s population has latent TB”, as one can find in the website of the World 
Health Organization (WHO) [49]. The value for the fraction of people infected with 
HIV is assumed ~ 1.7%, based on HIV & AIDS Information from AVERT.org [1]: 
“ There is either a generalised or concentrated epidemic. In a generalised epidemic, 
HIV prevalence is 1% or more in the general population. In a concentrated or low 
level epidemics, HIV prevalence is below 1% in the general population but exceeds 
5% in specific at-risk populations like injecting drug users or sex workers, or HIV 
prevalence is not recorded at a significant level in any group." The remaining values 
are estimated by assuming that we are in a “controlled” situation, without large 
percentages in the groups of highest risk such as A, At and Ch- Our aim is to 
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find the optimal combination of the fraction of individuals Ith that take correctly 
HIV and TB treatment (u][) or take only TB treatment (u^), in order to minimize 
the number of individuals with AIDS and TB diseases At- Different approaches 
were used to obtain and confirm the numerical results. One approach consisted in 
using IPOPT (short for “Interior Point OPTimizer”, a software library for large 
scale nonlinear optimization of continuous systems) [47] and the algebraic modeling 
language AMPL (acronym for “A Mathematical Programming Language”) [17]. A 
second approach was to use the PROPT Matlab Optimal Control Software [32]. 
For more details we refer the reader to [41, 42], where the same optimization ap¬ 
proaches are used. In Figure 2 we compare the extremal dynamics Cjj, I^ H and Ap 
associated to the extremal controls it* and it?] with the dynamics of the model (21) 
with Ui(t) = p and 112 (f) = Q, which coincide with model (3). In this simulations 
we consider /3i = 0.6, /?2 = 0.1 and the rest of the parameters take the values of 
Table 2, which corresponds to Rq = 4.91159 (I?i = 4.91159, R 2 = 1.07437). We 
assume that the weight constants take the same value W± = W 2 = 50. Observe that 
the number of individuals with AIDS and TB diseases At decreases significantly 
when the control measures it]], uj are implemented, see Figure 2 (c). On the other 
hand, the number of individuals that stays in the class Ch increases in opposition 
to the number of individuals that have both infections HIV and TB, see Figure 2 (a) 
and (b). During approximately 40 years the optimal combination of the fractions 



Figure 2. Dynamics Ch, Ith and At for cost functional (23), 
/?i = 0.6, fa = 0.1, W\ = W -2 = 50 and parameter values from 
Table 2. 


of individuals Ith that take HIV and TB treatments simultaneously and only TB 
treatment is around 0.5 and 0.46, respectively, see Figure 3 (a) and (b). In Figure 3 
(c) we observe that the extremal controls satisfy the restriction (22). At the end 
of 50 years, the number of TB and AIDS induced deaths reduces 5% when the con¬ 
trols u*, U 2 are applied, see Figure 4. Since the HIV treatments have higher costs 
than TB treatment, we can consider that the weight constant W\ associated to the 
control u\ takes greater values than IF 2 . In this case, the fraction of individuals 
that take TB and HIV treatment u[ decreases and the fraction of individuals that 
take only TB treatment increases, compared to the previous case W\ = W 2 = 50, 
but the associated extremal dynamics Ch, 1 ^ and Ap behave similarly to the ones 
in the case W\ = W 2 = 50, see Figure 5. The extremal controls in Figure 5 (a) and 
(b) are not intuitive, since the fraction of individuals that take both HIV and TB 
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(A) u* (B) u* 2 (C) u\ + u* 2 

Figure 3. Extremal controls u\ and u 2 for cost functional (23), 
Pi = 0.6, p 2 = 0.1, W\ = W 2 = 50 and parameter values from 
Table 2. 



Figure 4. Total population N for cost functional (23), Pi = 0.6, 
p 2 =0.1, Wi = W 2 = 50 and parameter values from Table 2. 





(a) u* (b) u 2 (c) u\ + 

Figure 5. Extremal controls u\ and u 2 for cost functional (23), 
Pi = 0.6, P 2 = 0.1, Wi = 500, W 2 = 50 and parameter values from 
Table 2. 


treatments is very low. If we assume that our aim is to minimize the cost functional 


Ji{ui{-),u 2 {-)) = f 
Jo 


A(t) 


Wi 

~Y 


M0 + 


W2 

2 


dt, (27) 


with T — 10 years and no disease induced deaths (dr = dr a = cZa = 0), that is, we 
wish to minimize the number of individuals that have only AIDS A and have both 
AIDS and TB diseases At, the extremal controls behave in a more intuitive way. 
Since we assume that there is no disease induced deaths we consider that is more 
adequate to consider T = 10 instead of T = 50 years. Moreover, is this case, the 
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total population is constant. We observe that the fraction of individuals that take 
both HIV and TB treatment u\ takes the maximum value for more than 7 years, 
and during this time the extremal control vanishes, see Figure 6 . In this case, we 



Figure 6. Extremal controls z and u 2 for cost functional J 1; 
with /?! = 0.6, /3 2 = 0.1, W\ = W 2 = 50 and parameter values from 
Table 2. 


compare the behavior of the dynamics A, At, Ith and Ch for the following cost 
functionals Ji, J 2 and J 3 , with W\ = W 2 = 50, where 


MM-)) = [ 

Jo 

MM-)) = [ 

Jo 


A(t) + A T (t) + — u\ (t) 


A(t) 


A r(t) + ^uKt) 


dt , 


dt, 


(28) 


(29) 


that is, when both controls u\ and U 2 are applied simultaneously, or are applied 
separately. The number of individuals with only AIDS, A , is lower for cost functional 


0.9 - 
0.8 - 

0.7 - - 

0.6 - 

| 0.5- 
O 

0.4 - 
0.3 - 
0.2 - 
0.1 - 

0 I-,-,-,-,- 

0 2 4 6 8 10 

time (in years) 



(A) U! 


(b) U 2 


Figure 7. Extremal controls when u\ and zz 2 are applied separately. 


Ji and with extremal controls given by u\ and W 2 in Figure 6 (a) and (b), see 
Figure 8 (a). However, this is not the best strategy for the reduction of the number 
of individuals with both AIDS and TB diseases, At- In this case, the best strategy 
is to apply only control zti where u± takes the values given in Figure 7, see Figure 8 
(b). The best strategy for the reduction of the total number of individuals with 
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only AIDS and with both AIDS and TB, A + At, during the first six years is the 
one associated to the cost functional J 2 , that is, apply only control u\ (see Figure 7 
(a)) and after six years the best strategy is to apply simultaneously both controls 
u\ and U 2 given by Figure 6. In Figure 9 we observe that the implementation of 



(a) A (b) A t (c) A + A t 


Figure 8. Dynamics A, At and A + At for cost functionals J\, 

J 2 , J 3 and u\ = p = 0.1, U 2 = q = 0.3, with = 0.6, ^2 = 0.1, 

W-[ = W 2 = 50 and parameter values from Table 2. 

controls ui and U 2 simultaneously or separately, contribute to the reduction of the 
number of individuals infected with HIV and TB, Ith, and increase the number of 
individuals that remain in the class Ch , that is, the HIV infection does not evolve 
to AIDS disease. In this case, the strategy of treating only TB for the individuals in 
the class Ith does not allow equal or better results on the reduction of individuals 
in the class Ith, which happens in the situation described in Figure 5 for the case 
W-! = 500 and W 2 = 50. 




Figure 9. Dynamics Ith and Ch for cost functionals J±, J 2 , J 3 
and ui = P = 0.1, 112 = q = 0.3, with /3i = 0.6, /?2 = 0.1, W\ = 
W 2 = 50 and parameter values from Table 2. 


7. Final comments and future work. Our numerical results only give extremals 
and no claims about optimality are made. As future work, it would be interesting 
to verify optimality by using second order optimality conditions and addressing 

















































A TB-HIV/AIDS COINFECTION MODEL 


23 


Symbol 

Description 

Value 

References 

N 

Total population 

variable 


N( 0) 

Initial population 

30000 


A 

Recruitment rate 

430 


A 4 

Natural death rate 

1/70 


Pi 

TB transmission rate 

variable 


02 

HIV transmission rate 

variable 


Vc 

Modification parameter 

0.9 


VA 

Modification parameter 

1.05 


ki 

Rate at which individuals leave Lt class by 
becoming infectious 

1/2 

[7, 21] 

Tl 

TB treatment rate for Lt individuals 

2 

[7, 21] 

T2 

TB treatment rate for It individuals 

1 

[7, 21] 

P'l 

Modification parameter 

0.9 


dr 

TB induced death rate 

1/10 

[7] 

5 

Modification parameter 

1.03 


v> 

Modification parameter 

1.07 


0 

HIV treatment rate for Ih individuals 

1 


pi 

Rate at which individuals leave Ih class to A 

0.1 


a l 

AIDS treatment rate 

0.33 

[3] 

Wl 

Rate at which individuals leave Ch class 

0.09 


d,A 

AIDS induced death rate 

0.3 


p2 

Rate at which individuals leave Ith class 

1 


V 

Fraction of Ith individuals that take HIV and 

TB treatment 

0.1 


<1 

Fraction of Ith individuals that take only TB 
treatment 

0.3 


T3 

Rate at which individuals leave Lth class 

2 



Rate at which individuals leave Lth class by 
becoming TB infectious 

1.3 fci 


r 

Fraction of Lth individuals that take HIV and 

TB treatment 

0.3 


P 2 

Modification parameter 

1.1 


CO 2 

Rate at which individuals leave Rh class 

0.15 


OL2 

HIV treatment rate for At individuals 

0.33 


dr a 

AIDS-TB induced death rate 

0.33 



Table 2. Parameters of the TB-HIV/AIDS model. 


properly the issue of conjugate points. For that, one needs to extend the theory 
underlying the computation of conjugate points and verification of optimality as 
developed in [38, Section 5.3] and [39]. 
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